function [mu, S2] = testGPReg(t)
%% GP Regression
%
% Tobias Siegfried, 10/01/09

figure(99)

xes = [-2:.2:2];
[xo,y] = meshgrid(xes);
xstarT = (-2:.2:2)'; xstar1 = repmat(xstarT,length(xstarT),1);
xstar2 = xstarT'; xstar2 = repmat(xstar2,length(xstarT),1);
xstar2 = reshape(xstar2,size(xstar2,1)*size(xstar2,2),1);
xstar = [xstar1 xstar2];
z = xo.*exp(-xo.^2-y.^2); %*randn(length(x),length(y));
figure(1)
subplot(2,2,1)
surfc(xo,y,z)

% pick some samples randomly
% threshold
% t = 0.5;
mRand = rand(length(xo),length(y));
mRand(mRand>t) = 1;
mRand(mRand~=1) = 0;
subplot(2,2,2)
surfc(xo,y,z.*mRand)
[x1,x2] = find(mRand);
x = [xes(x1)' xes(x2)'];
xi = sub2ind(size(mRand),x1,x2);
yRes = z(xi);
n = sum(mRand(:));
covfunc = {'covSum', {'covSEiso','covNoise'}};
loghyper = [log(1.0); log(1.0); log(0.1)];
subplot(2,2,3)
stem3(x(:,1), x(:,2), yRes,'+')
% compute interpolation
[mu S2] = gpr(loghyper, covfunc, x, yRes, xstar); % predictions
S2 = S2;% - exp(2*loghyper(3));
f = [mu+2*sqrt(S2);flipdim(mu-2*sqrt(S2),1)]; % standard noise-free pointwise errorbars
%fill([xstar; flipdim(xstar,1)], f, [7 7 7]/8, 'EdgeColor', [7 7 7]/8);
%hold on
%plot(xstar,mu,'b-','LineWidth',1);
%plot(x, yRes, 'k+', 'MarkerSize', 17);
%axis([-8 8 -3 3]);
subplot(2,2,4)
surfc(xo,y,reshape(mu,length(xo),length(xo))')
hold on
stem3(x(:,1), x(:,2), yRes,'k+')
hold off